Variational approximations to homoclinic snaking in continuous and discrete systems 
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Localised structures appear in a wide variety of systems, arising from a pinning mechanism due to 
the presence of a small-scale pattern or an imposed grid. When there is a separation of lengthscales, 
the width of the pinning region is exponentially small and beyond the reach of standard asymptotic 
methods. We show how this behaviour can be obtained using a variational method, for two systems. 
In the case of the quadratic-cubic Swift-Hohenberg equation, this gives results that are in agreement 
with recent work using exponential asymptotics. Secondly, the method is applied to a discrete system 
with cubic-quintic nonlinearity, giving results that agree well with numerical simulations. 
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I. INTRODUCTION 

This paper is concerned with the phenomenon of pin- 
ning of fronts in nonlinear dynamical systems. For sys- 
tems that exhibit bistability of two uniform states, a front 
connecting these two states will normally drift in one di- 
rection, depending on which of the two states is preferred. 
At a particular parameter value, known as the Maxwell 
point, there is no preference between the two states and a 
stationary front exists. However, in systems where there 
is an underlying structure, on a scale that is typically 
small compared with the lengthscale of the front, a front 
can become locked to this structure. This mechanism 
allows a stationary front to exist over a range of param- 
eter values around the Maxwell point, known as the pin- 
ning region. Placing two fronts back-to-back creates a 
localized state, and the bifurcation diagram plotting the 
length of this localized solution against the control pa- 
rameter has a snaking structure, involving a sequence of 
saddle-node bifurcations in near-perfect alignment. Since 
the spatial structure of such a localized state departs from 
and then returns to, a uniform state, the phenomenon has 
become known as "homoclinic snaking" . See P, [1] for a 
review of the subject and some open questions. 

There are two distinct scenarios in which homoclinic 
snaking appears. In the first of these, the underlying 
structure is provided by pattern formation. A physical 
system, for example convection in fluids 0, 01 , a vibrated 
granular material, buckling of a solid cylinder [5| , optical 
systems @ or gas discharge experiments 0, @|, has an 
instability leading to the formation of a regular periodic 
pattern as a parameter is varied. If this bifurcation is 
subcritical, there is bistability between the uniform and 
the patterned states. Near the onset of pattern forma- 
tion, a subcritical Ginzburg-Landau equation can be de- 
rived [1, 1^, [13] , which has front-like solutions connecting 
the two states. However, this equation does not capture 
the locking mechanism. The localised states have been 
found numerically in many studies [TT| - [l3j that use the 
Swift-Hohenberg equation as a simple model for pattern 
formation. 

A second situation in which locked fronts appear is in 
discretized forms of partial differential equations, where 
there is a locking effect to the imposed lattice. Examples 
include stationary solutions of the discrete bistable non- 



linear Schrodinger equation Il4l4l6l , which leads to a sub- 
critical AUen-Cahn equation |l7f, optical cavity solitons 
1^ Jj},], and discrete systems with a weakly broken pitch- 
fork bifurcation [23|. This situation also arises whenever 
a bistable system is solved numerically on a spatial grid. 

An interesting and challenging question is to determine 
the width in parameter space of the region in which sta- 
tionary localized solutions exist. Numerical simulations 
indicate that this pinning region becomes very small as 
the separation of the two length-scales increases [13, [13 ■ 
In fact the region is exponentially small, or beyond all 
orders. This is related to the fact that a conventional 
multiple-scales asymptotic method cannot describe the 
locking effect, since it regards to the two lengthscales 
as independent, while in fact the locking mechanism in- 
volves an explicit interaction between the short and long 
lengthscale. The necessary exponential asymptotics cal- 
culations, involving truncating a divergent asymptotic se- 
ries at an optimal point and obtaining an equation for 
the exponentially small remainder term, have recently 
been carried out for the Swift-Hohenberg equation with 
quadratic-cubic nonlinearities [2ll [22| and cubic-quintic 
nonlinearities [23]. However, these calculations are ex- 
tremely cumbersome and include an undetermined con- 
stant that must either be obtained numerically either by 
fitting to numerical results [2l| or by a ppr oximately solv- 
ing a complicated recurrence relation [23|. 

In this paper we use variational methods to obtain scal- 
ing laws for the structure of the snaking region, building 
on the work in a previous short paper that studied the 
cubic-quintic Swift-Hohenberg equation [23] . Of course, 
not all systems that exhibit homoclinic snaking have a 
variational structure, but most of those previously stud- 
ied do. The variational method is dependent on good ini- 
tial ansatz, but this is known in the cases studied here. 
For pattern forming problems near onset, it is known 
that the solution can be described by a slowly varying 
envelope function multiplied by a sinusoidal wave. Fur- 
thermore, the form of the envelope function is known 
from the standard asymptotic analysis, and can be ob- 
tained to higher order if required. Thus there is no guess- 
work involved in the method. Calculating the integrals 
in the Lagrangian automatically leads to exponentially 
small terms, and from these we can obtain the phase of 
the locked states that is inaccessible to the usual asymp- 
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totic expansion. Note that variational methods have been 
used before to study localized states on lattices [l^ , 
but not in the slowly- varying regime that is studied here. 
Also in the continuous case, variational methods have 
been used ^25j], but not in the snaking regime. 

Here, we consider two equations, namely the quadratic- 
cubic Swift-Hohenberg equation and the discrete cubic- 
quintic Schrodinger equation, also known as the spatially 
discrete AUen-Cahn equation, representing continuous 
and discrete systems, respectively. The first equation 
is discussed in Section |lll The results are then com- 
pared with numerical results obtained by continuation in 
Section Hi CI The second equation is studied in Section 
mil The scaling calculated analytically using the varia- 
tional method is then compared with computational re- 
sults, where good agreement is obtained. Conclusions are 
in Section IV. 



II. THE QUADRATIC-CUBIC 
SWIFT-HOHENBERG EQUATION 

The quadratic-cubic Swift-Hohenberg equation is given 

by 

dtu = ru^ {l + dlfu + b2u'^ -bsu^. (1) 

This equation represents a simple model for pattern 
forming systems that do not have a symmetry un- 
der sign reversal of the dependent variable u, and has 
been very widely used to illustrate homoclinic snaking 
Sinilillll. The Lagrangian for © is 

(2) 

It can easily be shown that £ decreases with time, so 
stable stationary states of ^ correspond to minima of 
(|2]). It will be assumed that 63 > 0; in this case u can 
be rescaled to set 63 = 1. In ([1]), the bifurcation at 
r = is subcritical (allowing localised patterns) if 62 > 
V27fe3/38. 

Figure[l]shows the bifurcation diagram for ([T]) obtained 
by a numerical continuation method for 62 = 1.5, 63 = 1, 
with periodic boundary conditions in a domain of length 
I = 827r. The norm N plotted is defined by 

.1/2 

^ dx. (3) 

J-l/2 

The periodic solution (dashed line) bifurcates subcriti- 
cally and becomes stable at a saddle-node bifurcation at 
r « —0.186. Two branches of localized solutions bifur- 
cate from the periodic state at a very small value of r and 



form an intertwined snaking pattern near the Maxwell 
point at r « —0.151. On one of these branches, the max- 
imum of the envelope function coincides with a minimum 
of the periodic pattern, and on the other it coincides with 
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FIG. 1: Bifurcation diagram of HI for 62 = 1-5, 63 = 1. The 
dashed line is the periodic solution. The two solid lines are 
the localized snaking solutions. 



a maximum, so the two states have a phase difference of 
TT. Solutions along the former branch are shown in Fig- 
ured The first figure, for r = —0.01, shows a slow spatial 
modulation that can be represented by an envelope in the 
form of a sech function (see Sec. Ill A| ). For r ~ —0.05 the 
amplitude modulation occurs over a shorter lengthscale 
but can still be represented by the sech function. The 
third graph shows u{x) at the first saddle-node bifurca- 
tion, r = —0.146, and the final graph is at a saddle-node 
bifurcation higher up the graph where r — —0.158 and 
N = 5.15. At this stage the solution resembles two fronts 
connecting the periodic solution to the state u = (see 
Sec.lTTB]). 



A. Analysis of solutions near r — 

The analysis near r — can be performed in a sim- 
ilar way to that in [l^l- Motivated by the results from 
multiple scale expansions [13, IH, , we take the ansatz 

u = A sech(i?a;) cos(fcx + ip) + C sech^ (Bx) , (4) 

from which we obtain that the norm ([3|) is given by 
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FIG. 2: Localized solutions on one branch of Fig.[T] From top 
to bottom, r = -0.01, r = -0.05, r = -0.146, r = -0.158 
(iV = 5.15). 



= _L (3^2 ^4^2) , 



, C(B^ + fc^) cos(i^) sech — + ^Bk cos(2tp) csch — 
is-' V V 2n I \ B 



Substituting the ansatz dU into the Lagrangian ([2]) yields the effective Lagrangian 
2C2 



L 



eff 



315B 



30B ^ 



(3663(72 - 5662C - 1685^ + + 105 - 105r) 

^-2463(72 + 2052C - 7B'^ - 30fc2B2 + lOS^ - 15k^ + 30fc2 + 15(r - 1)) + —A'^ 

SB 



(5) 



60b2C^k^B^ + 45 A% C B^k^ + B^ (240Cfc'* - ISOA^^a + 720C(1 - fc^)) + 2bs.C^k^^ 



7205^ ^ 

X cos(i^) ^e^*& + O ( e^'B' 



(6) 



It is immediately clear from (|6|) that the phase tp is de- 
termined by exponentially small terms, since the param- 
eter B is expected to be small, with 1/B representing 
the lengthscale of the modulation of the pattern. Fur- 
thermore, ^ shows that steady states (extrema of ie//) 
exist if is a multiple of tt [26|. There are two distinct 
states, — Q and ip = tt, both corresponding to even so- 
lutions but with the maximum of the envelope function 
sech(i3a;) coinciding with a maximum or minimum of the 
wave cos(A:a;), as shown by the numerical continuation 
method in the previous section. Note that this result does 
not depend on our choice of ansatz. For any slowly vary- 
ing, even modulation function f{Bx), the integrals in the 
Lagrangian involve powers of f{Bx) and its derivatives 
multiplied by cosikx) cos{ip) ^ yielding an exponentially 
small quantity multiplied by cos((p), and hence stationary 
states with = 0, tt. In the symmetric cubic-quintic case 
p4i] . only even powers appear in the integrals, leading 



to integrals only involving cos(2/cx) cos(2</?), and hence 
stationary solutions with p = Q, 7r/2, tt, 37r/2. 

Applying the Euler-Lagrange formulation to the effec- 
tive Lagrangian 



0, 



(7) 



gives us a system of nonlinear equations for a, with a = 
A, B, C, fc, p, that make Q an approximate solution of 

Neglecting the exponentially small terms in the effec- 
tive Lagrangian (jS]), the Euler-Lagrange equations can be 
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solved perturbatively about r = to yield 



B. Analysis of solutions near the Maxwell point 



A 

B 
C 
k 
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For values of 62 in the neighbourhood of 



_ 3&3 



-r + 0((-rf/2), 



-r + 0{-r), 



462 



4&2 - 363 



r + 0((-r)3/2), 



(8) 

(9) 
(10) 

(11) 



It is important to note that the presence of C 7^ 
plays an important role in the calculations above, un- 
like the case of the cubic-quintic equation [23] where it 
was possible to include only the first term in Q. Tak- 
ing C — would result in the leading order expression 
of A independent of b2, which is incorrect. This oc- 
curs because when C = the cubic term in the La- 
grangian, which corresponds to the crucial symmetry- 
breaking quadratic term in ([T]), does not contribute to 
the Lagrangian integral except through an exponentially 
small term. Note that C — 62^^/2, a result that is also 
easily be obtained at second order in the asymptotic anal- 
ysis of the problem. A similar Lagrangian method was 
used by Wadee and Bassom [2^, using more terms in 
the ansatz ^ but with (p = 0. The results above for A 
and C are almost the same as those obtained using mul- 
tiple scale expansions [H, HH] • A slight difference arises 
as we have for simplicity omitted a second-order term in 
sech^{Bx) cos{2kx 4- 2(^); including this term would give 
exact agreement between the Lagrangian method and the 
asymptotic analysis. According to (jS]) the transition from 
a sub- to supercritical bifurcation occurs at = 3/4, 

which is very close to the true value of 27/38 [2^. 

The equation (fTTj) indicates a slight decrease in the 
wavenumber, corresponding to a slight increase in the 
wavelength of the patterns. As noted by Wadee and 
Bassom ^25*1, this is simply a consequence of the linear 
dispersion relation; can be obtained directly from 
the linear terms in ([1]) when r < 0. Hence it is not sur- 
prising that the same wavenumber correction was found 
in the cubic-quintic case [23 |. 

When neglecting the exponentially small terms, the 
phase-shift ip at this order is arbitrary, which also agrees 
with the multiple scales result. However, taking into ac- 
count the equation d^pCeff = 0, in which all the terms 
are exponentially small, p has to be a multiple of tt, as 
mentioned above. 

The ansatz @ is only appropriate for small values of 
|r|, away from the Maxwell point (see Figs. [T] and [5]). In 
the snaking region near the Maxwell point, the envelope 
of the localized states resembles two connected fronts. 
The following section introduces a suitable ansatz for this 
regime. 



&20 = 



(12) 



the bifurcation is only slightly subcritical and the 
Maxwell point is within reach of weakly nonlinear anal- 
ysis. In this case the quadratic-cubic Swift-Hohenberg 
equation has a front solution, which is approximately 
given by (22| 



= A 



M- 



with 



Am 



tm 




-rM r 



-7" A/, 

'3 

"176166;^^'"^^°^ • 



(13) 

(14) 
(15) 



In analytically estimating the width of the snaking re- 
gion in the equation using variational approximations, 
we will employ a front solution similar to (jl3l) . Consider- 
ing this solution, one can note that the oscillation wave 
number fc of the front changes in space. In the limits 
X ±00, fc ^- 1 and (1 - ajif (27734)-!) « 1. There- 
fore, we will fix fc = 1 to simplify the analysis. 

It turns out that as in the case of the sech solutions 
considered in the previous section, if only the leading 
order front solution (|13l) is used as the ansatz, the leading 
terms in the Lagrangian do not depend on 62, giving 
qualitatively incorrect results. We therefore adopt an 
ansatz of the form 



Acos{x + (f)) &2^^ 

Tf^^fnfpry + 2(i-i-e'S(i^i--^))' 



(16) 



The second term here arises at second order in the weakly 
nonlinear expansion ^25.] , included for a similar reason as 
discussed in the previous section. This term is responsi- 
ble for the upward displacement of the periodic pattern 
that is apparent in Fig. [2] To simplify the calculation we 
have set the value of the coefficient of this term in advance 
rather than leaving it as a free parameter. The asymp- 
totic analysis includes another second-order term, pro- 
portional to cos 2 (a; -I- </?) [25^], which is omitted partly in 
the interests of simplicity and partly because it is smaller 
by a factor of 9 than the term we have included. 

Evaluation of the Lagrangian integral is considerably 
more complicated than in the cubic-quintic case [2^ . 
This is partly because of the two-term ansatz and partly 
because the dominant terms come from odd powers of 
the square root function, requiring integrals that can- 
not simply be computed using residues. These integrals 
can be evaluated in terms of multiplications of gamma 
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functions of a complex argument, which in turn can be 
expanded using Stirhng's formula when B is small, giv- 
ing the anticipated exponentially small terms. There are 
a large number of these exponentially small terms, and 
the calculation is complicated by the fact that terms in- 
volving higher powers of u, which one might expect to 
give a smaller contribution, in fact all give exponentially 
small terms of the same order. 

To illustrate some of the steps in the calculation, con- 
sider the term u^dx, which is the defined norm ^ 
and one of the terms in the Lagrangian density ^ . Using 
the ansatz, this can be written as 



u dx 



(1 + 62^2/2 + e^(--^)) 



2b2A^ cos(0) cos(a::) cos(20) cos(2a;) 



(l + eB(^-i)) 



3/2 



1 + f,B(^L+x) 



dxill) 



after using the double-angle formula and the fact that 
the envelope function is even. The first term, that is in- 
dependent of 4>, in the integral can be evaluated directly 
and gives the answer A^L (l + ^^62/2) - A-^by{2B) + 
0{e~^'"). The integral of the third term can be found 
by using contour integration, similarly to that used in 
[2J|. Yet, we are not interested in the explicit expression 
of it, as it contributes to higher harmonic corrections in 
(j). The second term is of our interest, giving a leading 
exponentially small contribution of order 0{e~^). Un- 
fortunately, the integral cannot be immediately evaluated 
by using the residue method as mentioned above. By us- 
ing Mathematica, one will obtain that 



Jo 



cos(a;) dx 

[I + e.Bi^x-L) 



~iLj 



'■iB+2i 
\ 2B 



-i+B\ 



3/2 



^BL 



B 



a-BL 



aBL 



2F1 



B 1 



-BL 



B '2'B' 



B 



B 1 
~ 2' 



I 

B' 



-BL 



(18) 



where F and 2^1 are respectively the Gamma and the 
hypergeometric function. 

The hypergeometric function can be written in a series 
form as 



{a)nib)n 
(c)„ n! ' 



where (•)„ = •(• -I- 1) . . . (• -f n — 1) is the Pochhammer 
symbol. The series is convergent in our case, i.e. < 
B < 1 and Si > 1, such that up to ©(e^^SLj 



^F^ 



2F1 



i~B 1 i 



-,-BL 



1 



i + B 1 -i 



-,-BL 



B T B 



I - -e-^^ [l-iB). 



when z is large enough in absolute value. Then, the first 
two terms on the right hand side of (|18l) can be approxi- 
mated by 



((-8 + m cos(L) + (8 + 35) sin(L)). (19) 

Combining the approximations above, one will obtain 
that 



u'dx = A^L{l + -A%i\ - ^ 
' 2 2; 2B 



00 



X (-8 + 3B) cos(L) + (8 + 35) sin(L) . (20) 



As for the Gamma function, it can be approximated 
by the Stirling's formula 



r(^) 



z \eJ \ z 



The formula (|20l) indicates that in the snaking region, so- 
lutions with a larger norm correspond to longer plateaus. 

Performing the same calculations to all the terms in 
the Lagrangian ^ yields the effective Lagrangian 



6 



eff 



-(15B'' + 480^2 



l4. .2 



1920B 

+2iOA% - 360A^b3 - 1080^^636^ - 110A%3b^) 
LA^ 



480^1^5^ 



ASOA^rb"^ 



96 



■(18^2^3 



2AA^bl 



48r - 2AA''rb^ - 8A%^ + 36A%3bi + 3A%3b^) 



' cos((^) 



420 



(e-5^(sin(2L),cos(2L)),e' 



BL 



V 



(21) 



where 



X± ±^^A%\b:i^2%^A?b\B 



±156^2 ±840 rB^ 



A20 A%B 
763B^ + 70 rB^. (22) 



Consider first the terms that are not exponentially 
small, that is, the terms that do not involve the phase 
(p. Bearing in mind the anticipated scaling from (1141) . 
A = ©(r^/*), the dominant terms in the bracket multi- 
plied by L are the first two terms, so Ceff only has a 
minimum over L if at leading order, 63 = 363/4. This 
is very close to the asymptotic result (fT2)) . and would be 
exactly the same if we had included the cos2(x + (p) term 
mentioned above. We therefore set 



62 = 7363/4 + A, |A|«1, 



(23) 



and expect that A and B are ©(r^/^), in which case the 
leading terms in the Lagrangian are 



nlead 



A^ 



USB 
LA^ 



{32B^ + 32^/3hA^A - 45A'*6^) 



+ -^(45^'*6^/2 - 48r - 24V363yl^ A)(24) 

Finding and solving the Euler-Lagrange equations by dif- 
ferentiating with respect to L, A and B gives the values 



A' = 



8\/3A 



156; 



3/2 ' 



B = 



lOA 



5Vh 



and yields the Maxwell point 

2A2 



rM = 



563 



= -B' 



(25) 



(26) 



which is within 3% of the asymptotic value given in (|15p . 
Note that the length L is not determined by these equa- 
tions. Substituting these values of A, B and tm into the 
leading terms of the Lagrangian we find that the mini- 
mum is 



4V30A2 
^56r 



(27) 



Thus, although the Maxwell point can be determined by 
the condition that £ is zero for both the zero state and the 



periodic state, the value of L is not zero for the localized 
state, due to the presence of the fronts. In fact, £ is 
positive and of the same order as r. 

Considering now the exponentially small term involv- 
ing (/?, it is clear from d^Ceff = that there are two 
distinct branches of stationary snaking solutions, with 
(/? = and </? = TT, corresponding to the centre of the 
localised pattern being a local maximum or minimum 
respectively. However, another possible way to satisfy 
d^Ceff = is that the coefficient multiplying the cost/? 
term vanishes, i.e. 



K+ cos(L) + A'„ sin(L) = 0. 



(28) 



Substituting the above values of A and i3, and taking 
the leading terms in K± , solutions of this type occur for 
values of L given by 

timL = K+/K^, (29) 

where K± is given in \22\ . which to the leading order is 



8^0F 
7875 63 



(±307 VS- 210 



(30) 



Hence, for small B these states exist if L = — 0.110219. ..+ 
TOTT for integer values of m. These solutions have been 
referred to as "bridges" [271 or "ladders" |23| and were 
identified in the exponential asymptotics of ([T]) by Yang 
and Akylas [2^ and Chapman and Kozyreff |21|]. The 
ladder states (not shown in Fig. [T]) have no reflection 
symmetry and connect the two snaking branches, bifur- 
cating from them near the saddle-node bifurcations. 

To find the snaking range, we set sin((^) = and r = 
tm + and the Lagrangian simplifies to 



L 



4V30A2 6rA^L 



7561 



+e~* (ii:+cos(L) +ii:_sin(L)), (31) 



with K± given in pop . This can be minimised over L 
to give the exponentially small value of 5r for stationary 
solutions. 



5r 



2e- 



A^ 



{K_ cos(i) - K+ sin(X)) . (32) 
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with 5r = and (b) with 5r = 0.4. 

The maximum value of 5r (half the width of the snaking 
region) is then, using the leading-order approximations 
for 62, A and B given above, 

Sr,n = %^ 7241249^^. (33) 

Note that the dependence of the snaking width on the 
small parameter A is very similar to that for the small 
parameter 63 in the cubic-quintic case [l^- The depen- 
dence in (|33l) is the same as that obtained by Kozyreff and 
Chapman 22] using the methods of exponential asymp- 
totics, although the comparison is not immediate since 
they regard 62 as the control parameter rather than r. 
Note that the numerical coefficient was determined by 
fitting in 22]. 

A useful feature of the Lagrangian method is that we 
can examine the stability of steady states since we know 
that stable equilibria are local minima of the Lagrangian. 
Considering the exponentially small terms, with A and 
B fixed and r = tm + the dependence of £ on 5r, L 
and If) can be represented by the function 

Ce{5r, L, ip) — —SrL + cosipcos{L — Lq), (34) 



where the phase shift Lq satisfies (|29p and for convenience 
we have set the constants to 1. Figure [31[a) shows a sur- 
face plot of £e for 6r — 0. The solutions on the snaking 
branches, with ip — 0,?:, L — Lq — mn are either maxima 
or minima, so are either stable or unstable with two pos- 
itive eigenvalues. The ladder solutions at ip = 7t/2, 3tt/2, 
L — Lq = 7r/2 + TOTT are saddle points and are therefore 
always unstable. The surface plot for Sr = 0.4 is shown 
in Fig. [SJb) . The symmetric snaking solutions remain at 
= 0, TT but are shifted in L, while the ladder states are 
shifted in (p but remain at the same values of L. 

C. Numerical results 

In this section we compare the above results from 
the variational method with numerical solutions of the 
quadratic-cubic Swift- Hohenberg equation ([T]). We have 
solved the equation numerically for localized states, us- 
ing a pseudo-arclength continuation method with peri- 
odic boundary conditions, implemented with a Fourier 
spectral discretization. Plotted in the top left panel of 
Fig. m is the bifurcation diagram showing two branches 
of localized solutions for &2 = 1-3 and 63 = 1. In the 
same panel, shown in dashed lines are our analytical re- 
sults obtained from ([T]) using the sech ansatz, where one 
can see that the variational calculation approximates the 
numerics very well for relatively small |r|. As r decreases 
further and enters the snaking region, the approximation 
deviates from the numerics. In the top right panel, we 
plot the profile of a localized solution at r = —0.06 and 
its approximation. 

Considering Eqs. (j8|)~ ([Tl]) one can conclude that for 
larger values of the parameter 62, the parameter |r| can 
be larger while keeping the amplitude A ^ small. This 
implies that the sech ansatz can have a longer validity 
region for large 62 ■ It is important to note that the ansatz 
might be able to quantitatively capture the first snaking 
bifurcation for a larger value of 62. 

In panel (c) of the same figure, we plot the Lagrangian 
^ corresponding to the bifurcation diagram (a) as a 
function of the length of the solution's plateau 2L, which 
is calculated numerically as 



{UM - Urn)"^ -I- 2{UM + Um)^ ' 

where um = max{u} and Um = niin {u}. The integration 
is a definite integration over the computational domain. 
Plotted in the same panel is our effective Lagrangian pip . 
There is good agreement for the average numerical value 
of the Lagrangian and the qualitative nature of the os- 
cillations, but the amplitude of the oscillations appears 
to be underestimated in the variational approximation. 
Note that the amplitude of the oscillations in C increases 
with L since from ((3T|) with 5r oc cos(L), there are oscil- 
lations of the form Lcos(L). 

Finally, in panel (d) we show the width of the snaking 
region as a function of 62 numerically and analytically. 
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FIG. 4: (a) The bifurcation diagram obtained numerically and our approximation obtained from solving ^ with the ansatz (Q 
for 62 ~ 1-3. (b) A numerically obtained localized solution and its approximation from the variational formulation for (p — 
and r — —0.06. (c) The numerically computed Lagrangian ((2| as a function of the length of the plateau 2L, corresponding to 
panel (a). The dashed line is our approximation calculated from the effective Lagrangian (|31|l . (d) The width of the snaking 
region as a function of 62. Filled circles are numerical and solid lines are analytical, i.e. 2Srm (|33p . The inset shows the same 
comparison in a log scale. In all the figures, 63 = 1. 



where one can see that our approximation is in good 
agreement. 

The numerical method can also be used to check the 
stability of the steady state solutions (since the Jacobian 
is already computed as part of the Newton-Raphson it- 
eration). This confirms that the snaking states are in 
general either stable, or unstable with two small positive 
eigenvalues. One eigenvalue changes sign at the saddle- 
node bifurcation and another at the nearby bifurcation 
to the ladder states. The ladder states are found to be 
unstable with one positive eigenvalue, consistent with the 
fact that they are saddle points of the Lagrangian. 



III. SNAKING IN DISCRETE SYSTEMS 

To illustrate the wide-ranging applicability of the vari- 
ational method in describing the snaking behaviour, we 
analyse in this section an analogous problem in discrete 
systems. This turns out to be remarkably similar to the 



continuous case. Consider the system of ordinary differ- 
ential equations 



dun 
IT 



C{u 



n+l 



2Un + Un^l) + ^lUn + 2ul^ - ul, (35) 



with the condition that u„ — >■ as n — >■ ±00. We assume 
that C > and look for stationary solutions of ([55]) . 
This system has been considered in l4-17 , and very 
similar systems have been studied in 1^ |20l • The sys- 
tem ([55)1 is interpreted as a discrete form of the nonlinear 
Schrodinger equation in (Tsl 16 [ , and a discrete subcriti- 
cal AUen-Cahn equation in [17| . These previous studies 



have used numerical continuation to obtain bifurcation 
diagrams that show a very similar snaking structure to 
that seen in the continuous Swift-Hohenberg equation 
((H). Note however that ([55t is not a discrete form of 
the Swift-Hohenberg equation. Rather it is a discrete 
form of the subcritical Ginzburg-Landau equation that 
describes the slowly varying envelope function for ([I}. 
Thus the snaking behaviour seen in (j35p results from a 
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locking mechanism to the discrete lattice, with the lattice 
points providing the analogue of the small-scale pattern 
in the Swift-Hohenberg equation. In view of this, com- 
putational investigation of the snaking in (1351) is much 
easier than in ([ij, since there is no small-scale structure 
to resolve. 

The grid-locking behaviour of fronts in (j35|) can be seen 
as an indication of numerical error in the finite-difference 
approximation of the bistable partial differential equation 



du 
'dt 



(36) 



When is discretized with second-order finite dif- 

ferences with mesh spacing e, the result is (1351) with 



C — e~^. Front-like solutions to (1551) are travelling waves, 
for all values of /i except exactly at the Maxwell point. 
Thus the snaking branch of (I55t indicates a range of val- 
ues of where the numerical approximation incorrectly 
finds a stationary solution. It is therefore of interest to 
determine the width of this region for small values of e 
(large values of C). An alternative scaling, similar to that 
in Sec. [Ill is to set C = 1 and replace the 2 in (I55|) by a 
small parameter 62 ; however we will use the form (1351) to 
facilitate comparison with the previous work cited above. 
The Lagrangian for (l35t is 



C 



M 2 
2"" 



1 



6 



(37) 



Note that a variational method has been used in jlj, |15| 
to study the case of small C, but here the focus will be 
on large C . The uniform solutions of p5p are given by 
Un = and 



1± 



(38) 



so there is a saddle-node bifurcation at ^ = — 1 and there 
is bi-stability of both the zero state and the larger non- 
zero solution in (P5| for — 1 < /i < 0. By finding the 
Lagrangian of this non-zero state it is straightforward to 
show that the Maxwell point is at /i = —3/4 and so at 
this point the stable non-zero state is it„ = ■\/3/2. 

For large C we expect the solution to ((35)) to be close 
to that of the continuous system ([5S)) . At the Maxwell 
point there is an exact stationary front-like solution of 
\, u{X), given by 



1 



, where = 3/2, a = Vs. 



(39) 



In the discretized system, m„ = u{en), and we wish to 
consider a localized state linking two fronts. The ansatz 
analogous to that used in the continuous case (jl6p is then 



I _(_ gae(|n-0|-L) ' 



(40) 



Here, the length of the front (measured in terms of the 
number of mesh points) is 2L and is a phase variable, 



which is arbitrary at this stage but will be determined 
by exponentially small terms. If = then Uj = U-j 
and the localized state is symmetrical and 'site-centred'. 
Similarly, a 'bond-centred' symmetrical state, with the 
property Uj — ui-j, can be obtained by setting (j) = 1/2. 
Note that it is not necessary to include higher order terms 
in (HD|) as in ([T5|) ; this is because ([55]) only has nonlinear 
terms with odd powers and therefore is more analogous 
to the cubic-quintic Swift-Hohenberg equation analysed 
with the Lagrangian approach in [24] . 

Consider now the sum J2'^=-oo tiia.t appears in C. 
This can be evaluated by writing the sum as an integral 
of a function multiplied by a sum of 6 functions, and then 
writing the sum of 6 functions as a Fourier series: 



A^ 



— — 00 

A^ 



I _|_ gae{\x-4>\-L) 
A^ 

I -|_ gae{\x~(t>\~L) 



S{x — n) dx 



n— — 00 
00 

E 

k——oo 



dx. 



Note that x here is not the same as X in (pG]) . rather it is a 
continuous form of the variable n. Of course the Fourier 
series representation of the 6 function is not uniformly 
convergent, but the integral converges very rapidly - in 
fact exponentially. From the term k — we obtain an 
order one contribution 



A^ 



1 _(_ gQe(|a:-0|--L) 



dx = 2LA^ + 0(e-"'-^). (41) 



The terms with fc = ±1 give a contribution 



1 -)_ Qae{\x-ct>\-L) 



dx 



(42) 



which is exponentially small in e, and the terms arising 
from larger values of k are smaller by an exponentially 
small factor. Note that the integrals that are needed here, 
for the discrete case, are essentially exactly the same in- 
tegrals that were needed for the continuous case in Sec.HIl 
and 0. 

The integral (|42|) can be found using contour inte- 
gration. For the term in e^"^ the contour is closed in 
the upper half plane and is dominated by the poles at 
X = (f) zL L + in/ (ae). The sum of the residues at these 
two poles is 



ae 



(43) 



The value of (j42p is found by multiplying this by 2i7r and 
adding the complex conjugate to account for the e^'^*'^^ 
term. Hence the required sum is 
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ul = 2LA^ + 8Trsm{27rL)cos{2TT(j>) — + 0(e^ , e""'-^). (44) 



ae 

n— — oo 



As in the continuous snaking case, we suppose that the exponential terms responsible for the grid locking dominate 
those from the interaction between the two fronts. This requires L e^^. 
In a similar way it can be shown that 

°° 24^ HnA'^ -2 2 

> ut^2LA^ ^cos(27r</>)(27rcos(27rL) - aesin(27rL))e^^ (45) 



ae a^e 

n— — oo 



and 



V u*5 = 2LA'^ 5-^cos(27r<^)(37raecos(27ri) + (2n^ - a^e^) sin(27rL))e— , (46) 



n— — oo 



where again terms of order e and e""*^^ have been dropped. Note that for small e, the dominant exponentially 
small term comes from (j46p . 

The remaining term required for the Lagrangian is 

oo oo ^ \ 2 

^ {Un+l - Unf ^ Y L ^gae(|n+l-0|-L) " 1 + e^^d"--^!-^) j ' ^'^''^ 

n— — oo n— — oo ^ ^ 

As with the other terms, this can be written as an integral, 

^ f 1 1 \ ^ 

\l+e°'<\^+^-'f'\-L) l_^gaei\x^cl>\-L) J ^^^^ +^ + . . .j OX, 



leading again to a part that is independent of and an ex- 
ponentially small 0-dependent term. For the part which 
does not depend on </), we can write the inte grand as 
(/(X + e)-/(X))2, where /(X) = A/Vl + e"^, X = ex, 
and use Taylor expansion to find the dominant contribu- 
tion to (H71) in the form 

\A^ae-^^AW + Oie% 

If the same method is applied for the exponentially small 
terms, it turns out that each term in the Taylor expansion 
yields a contribution that is of the same order after the 
exponentially small integral has been evaluated. Each 
term gives a contribution proportional to 

— cos(27r(/)) sin(27rL)e , 
ae 

with a different numerical constant. Thus we know the 
scaling but not the magnitude of this term. Since (|47p is 
multiplied by C = in (l37t . this term is of the same 
order as the term in ([571) . 

Using the above results, the leading terms in the La- 
grangian ([571) are 

A^a A^ / . . A^\ 

8e ae 2ae \ A j 



By making this expression stationary with respect to A, 
L and a we recover the results given earlier, A^ — 3/2, 
a — -v/S, = —3/4. 

To find the width of the snaking region, we set fi — 
— 3/44-A^, and introduce the largest of the exponentially 
small terms, arising from the term and the difference 
term in (|37l) . to obtain 

C=^ T^LAfi- "- / cos 27r(/. sin 2^£ e— , 

8e 2 oa'^e'^ 

(49) 

where D is an unknown constant representing the con- 
tribution from (|47p where only the scaling is known. 

By making ([1^ stationary with respect to (p it fol- 
lows that snaking branches branches exist with (f) = or 
(f) = 1/2, corresponding to the 'site-centred' and 'bond- 
centred' states found numerically. For these snaking 
branches, differentiating with respect to L shows that 

1 -2„^ 

Afi oc — cos(27rL)e 

so the scaling for the width of the snaking region in ([35]) 
for large C is 

cx c3/2e-2-'V^. (50) 
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FIG. 5; Log-log plot showing numerically obtained values of 
^^g2T ^Jcji (^pQjji^^g) together with the 3/2 scaling law ex- 
pected according to (|50p (line). 

Furthermore, differentiating (|49p with respect to pre- 
dicts that the ladder states, with sin(27r(/)) 7^ 0, occur for 
integer values of L. 

To check the scahng predicted by (j50p . snaking curves 
for the equation were obtained numerically using a con- 
tinuation software method. With 200 grid points it was 
possible to obtain snaking bifurcations up to a value 
C = 14 (at which point the snaking width is of the 
order of 10"^''). Figure [5] shows the numerically ob- 
tained snaking width multiplied by the predicted expo- 
nential factor e^'^ VcV3_ This shows a clear power-law 
behaviour, with a power very close to the value 3/2 pre- 
dicted by (Uni). 

IV. CONCLUSION 

In this paper we have used the variational approxi- 
mation to study the snaking behaviour of localised pat- 



terns in the quadratic-cubic Swift-Hohenberg equation 
and the discrete bistable AUen-Cahn equation. With a 
simple ansatz, inspired by asymptotic analysis, the expo- 
nentially small terms responsible for the snaking appear 
in the Lagrangian. This enables the branches of snaking 
solutions to be found, along with the asymmetric 'ladder' 
states that link these branches. 

These solutions cannot be found by a conventional 
multiple-scales method, since they involve a locking 
mechanism between the long and short scales, but are 
accessible through exponential asymptotics 21, 22]. The 
Lagrangian approach provides a useful complement to 
the exponential asymptotics method. Both methods give 
the same scaling for the relationship between the width 
of the snaking region and the small parameter of the sys- 
tem. 

We have shown that a close similarity exist between 
the pinning phenomena in continuous and discrete sys- 
tems. This arises partly through the fact that the con- 
tinuous limit of the discrete system considered is exactly 
the same as the Ginzburg-Landau equation describing 
spatial modulation of the pattern in the continuous case, 
and partly because the sums in the Lagrangian for the 
discrete case can be converted into integrals very similar 
to those that appear in the continuous problem. 

The variational method has a number of advantageous 
features, in addition to the fact that the Lagrangian 
integral immediately generates the necessary exponen- 
tially small terms. Based on the minimisation of the 
Lagrangian, it is easy to distinguish between stable equi- 
libria and unstable ones. Furthermore, although we have 
concentrated on the case of small parameters, for the 
purposes of comparison with asymptotic methods, the 
method is not restricted to this regime. In future work it 
may be possible to apply the method to cases where no 
parameters are small, and make useful comparisons with 
numerical or experimental results. An additional chal- 
lenge will be to extend the results to two-dimensional 
systems [Ullllil. 
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